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Integration of solar energy into the electricity network is becoming essential because of its continually 
increasing growth in usage. An efficient use of the fluctuating energy output of photovoltaic (PV) systems 
requires reliable forecast information. In fact, this integration can offer a better quality of service if the 
solar irradiance variation can be predicted with great accuracy. 

This paper presents an in-depth review of the current methods used to forecast solar irradiance in 


Keywords: 

Solar irradiance 
Forecast models 
Statistical models 

NWP models 
Postprocessing methods 


order to facilitate selection of the appropriate forecast method according to needs. The study starts 
with a presentation of statistical approaches and techniques based on cloud images. Next numerical 
weather prediction or NWP models are detailed before discussing hybrid models. Finally, we give 
indications for future solar irradiance forecasting approaches dedicated to the management of small- 
scale insular grids. 

© 2013 Elsevier Ltd. All rights reserved. 
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1. Introduction 


The contribution of photovoltaic systems (PV system) power 
production to the electric power supply is constantly increasing. Utility 
companies and transmission system operators have to deal with the 
fluctuating input from PV system energy sources. This is a new 
challenge compared with power production from conventional power 
plants that can be adjusted to the expected load profiles. An efficient 
use of the fluctuating energy output of PV systems requires reliable 
forecast information. 

Note that the reliable forecasting of the expected solar resource 
is but one aspect of the broad question of solar resource assess- 
ment that ranges from, for example the work of Perez et al. [1] on 
variability to Lucia [2] on the link between the entropy generation 
maximum principle and the exergy analysis of engineering and 
natural systems. In this paper, we focus solely on the, once again, 
broad field of forecasting, broad both in approaches taken, and the 
time scales covered. 

Load patterns forecasted for the next 2 days provide the 
basis for scheduling of power plants and planning transactions 
in the electricity market in order to balance the supply and 
demand of energy and to assure reliable grid operation [3]. These 
forecasts are used by utility companies, transmission system 
operators, energy service providers, energy traders, and indepen- 
dent power producers in their scheduling, dispatching and regula- 
tion of power. 

In particular, insular territories experience an unstable elec- 
tricity network and use expensive means in order to provide the 
power for the peak demand periods. Their grids are generally 
not interconnected with any continent and all the electricity 
must be produced inside the territory. The power of grid 
connected PV plants increases fast and can interfere with net- 
work stability. An efficient forecasting method will help the grid 
operators to better manage the electrical balance between 
demand and power generation. Kostylev and Pavlovski [4] 
identify three forecasting horizons (intra-hour, intra-day and 
day ahead) related to the grid operator activities (ramping 
events, variability related to operations, unit commitment, 
transmission scheduling, day ahead markets, hedging, planning 
and asset optimization). 

Forecasting of global horizontal irradiance (GHI) is the first and 
most essential step in most PV power prediction systems. GHI 
forecasting approaches may be categorized according to the input 
data used which also determine the forecast horizon. 


e Statistical models based on online irradiance measurements are 
applied for the very short term timescale from 5 min up to 6h 
(see Reikard [5]). Examples of direct time series models are 
autoregressive (AR) and autoregressive moving average 
(ARMA) models. Furthermore, artificial neural networks (ANNs) 
may be applied to derive irradiance forecasts. 

e For short-term irradiance forecasting, information on the 
temporal development of clouds, which largely determine 
surface solar irradiance, may be used as a basis. 

© Forecasts based on cloud motion vectors from satellite 
images (Lorenz et al. [6]) show good performance for the 
temporal range from 30 min up to 6h. 

o For the subhour range, cloud information from ground- 
based sky images may be used to derive irradiance forecasts 
with much higher spatial and temporal resolution compared 
with the satellite-based forecasts. 

e For longer forecast horizons, from about 4 to 6h onward, 
forecasts based on numerical weather prediction (NWP) mod- 
els typically outperform the satellite-based forecasts (see Perez 
et al. [1], Heinemann et al. [7]). 

e There are also combined approaches that integrate different 
kinds of input data to derive an optimized forecast depending 
on the forecast horizon. 


Solar irradiance forecasts was assessed in terms of root mean 
square error (RMSE) and mean bias error (MBE or bias) which are 
defined as follows: 


1 


i 


n 
RMSE = bs (Xpred.i—Xobs.i)” a) 
=1 


a 


MBE = 2 x (Xpred,i—Xobs i) (2) 
i= 

where Xpreqj ANd Xops; represent the ith valid forecast and observa- 

tion pair, respectively and n is the number of evaluated data pairs. 

These metrics are not formulated in the same way in all the papers 

we reviewed. David et al. [8] illustrated several formulas wrongly 

called RMSE or MBE. 

Many solar irradiance forecasting models have been developed. 
These models can be divided into two main groups: statistical 
models and NWP models. Statistical models are based upon the 
analysis of historical data. They include time series models, satellite 
data based models, sky images based models, ANN models, wavelet 
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analysis based models, etc. NWP models are based on the repro- 
duction of physical phenomenon. 

The paper is organized as follow. In Section 2, statistical 
approaches are presented. In Section 3, cloud imagery and satellite 
based models proposed in the literature are reviewed. In Section 4, 
the NWP approaches presented in the literature are reviewed. In 
Section 5, hybrid models are evaluated. Finally Section 6 is 
dedicated to trends for future solar irradiance forecasting in an 
insular environment. 


2. Statistical models 


Forecasting methods based on historical data of solar irradiance 
are two categories: statistical and learning methods. Seasonality 
analysis, Box—Jenkins or Auto Regressive Integrated Moving Aver- 
age (ARIMA), Multiple Regressions and Exponential Smoothing are 
examples of statistical methods, whilst AI paradigms include fuzzy 
inference systems, genetic algorithm, neural networks, machine 
learning, etc. 


2.1. Linear models or time series models 


Statistical methods have been used successfully in time series 
forecasting for several decades. Using the statistical approach, 
relations between predictors, variables used as an input to the 
statistical model, and the variable to be predicted, are derived 
from statistical analysis. Several studies with respect to direct time 
series modeling have been performed. In Reikard [5], different 
time series models are compared. In Bacher et al. [9], the authors 
investigate the use of a simpler AR model to directly predict PV 
power in comparison with other models. 


2.1.1. Persistence forecast 

It is useful to check whether the forecast model provides better 
results than any trivial reference model. It is worthwhile to imple- 
ment and run a complex forecasting tool only if it is able to clearly 
outperform trivial models. Probably the most common reference 
model in the solar or wind forecasting community for short term 
forecasting is the persistence model. The persistence model sup- 
poses that global irradiance at time t+1 is best predicted by its value 
at time t 


Xen =X; 


The persistence forecast, also known as the naive predictor, can be 
used to benchmark other methods. Persistence forecast accuracy 
decreases strongly with forecast duration as cloudiness changes from 
the current state. Generally, persistence is an inaccurate method for 
more than 1h ahead forecasting and should be used only as a 
baseline forecast for comparison to more advanced techniques. 

In Perez et al. [1], the single site performance of the forecast 
models is evaluated by comparing it to persistence. 


2.1.2. Preprocessing of input data 

When using statistical time series analysis, any type of conditional 
forecast model is structured to deal with stationary series, at least 
weakly stationary. This means no trend nor seasonality, and the 
series is homoscedastic (constant variance). There are several ways to 
deal with non-stationary series to get them into an appropriate form. 

Processes to obtain stationary solar irradiance time series: The 
solar insolation is the actual amount of solar radiation incident 
upon a unit horizontal surface over a specified period of time for a 
given locality. It depends strongly on the solar zenith angle. For 
statistical models, it may be favorable to treat the influences 
of the deterministic solar geometry and the non-deterministic 
atmospheric extinction separately. For this purpose, two 


transmissivity measures have been introduced: clearness index 
(k) and clear-sky index (k*). 

Clearness index. The clearness index k is defined as the ratio of 
irradiance at ground level I to extraterrestrial irradiance Text on an 
horizontal plane 


k=I/lexe (3) 


It describes the overall extinction by clouds and atmospheric 
constituents in relation to the extraterrestrial irradiance. This 
approach strongly reduces seasonal and daily patterns by consider- 
ing the influence of the zenith angle, which is modeled by I,,;. The 
clearness index is widely applied to reduce the deterministic trend 
in irradiance time series. However, the clearness index accounts for 
only the trends caused by geometric effects on solar position. As 
atmospheric extinction depends on the length of the path of the 
radiation through the atmosphere, it is also governed by solar 
geometry [3]. 

Clear sky index. The clear-sky index decreases with increasing of 
solar zenith angle. To account for this influence as well, the clear- 
sky index k* is introduced. k* is defined as the ratio irradiance at 
ground level I to irradiance of a defined clear-sky model I¢jear 


k* = I/Tctear (4) 


For the calculation of the clear-sky index, a clear-sky model and 
information on atmospheric input parameters are required. A clear 
sky model estimates the global irradiance, usually referred to as clear- 
sky irradiance Iciear, in clear sky conditions at any given time. An 
overview of different models is presented by Ineichen [10]. Clear-sky 
models range from empirical models to radiative transfer-based 
calculations. All these models need information on the state of the 
atmosphere as input [3]. 

The quantities introduced in this section are frequently used in 
solar modeling and forecasting. For example, some time series 
models explicitly require input parameters free of trend; hence, 
clearness index or clear-sky index can be an adequate choice. Also, 
satellite-based forecasts of irradiance are based on the concept of 
separately describing the influence of clouds and other atmospheric 
components by using the clear-sky index and a clear-sky model. 
Furthermore, most empirical models to derive the diffuse fraction of 
irradiance, necessary to calculate the irradiance on a tilted plane, are 
generally based on the clearness or clear-sky index. 

Preprocessing of the input data can considerably contribute to 
improving the accuracy of forecasts, and different approaches are 
proposed. As mentioned earlier, stationary, trend-free time series are 
required for classical time series approaches, and might be beneficial 
also for ANN. Hence, the use of the clear-sky or clearness index 
instead of irradiance data seems suitable. This approach is followed, 
for example, in Bacher et al. [9], and Kemmoku et al. [11]. On the 
other hand, Sfetsos and Coonick [12] argues that time series of the 
clearness or clear-sky index are mostly random, and hence do not 
provide a good basis for any learning algorithm. They recommend 
using irradiance values as input. Other examples of preprocessing of 
input data are the use of wavelets in Cao [13] and the use of the 
logarithm of irradiance values in Reikard [5]. 

Statistical tools for removing trend and seasonality. Because of its 
unpredictable noise, it is not easy to find the trend in a day's series 
of solar irradiance. Several models exist to detrend the hourly solar 
irradiance (Baig et al. [14], Kaplanis [15]). Clear sky index, as 
described above, is one model to obtain the deterministic daily 
variation of irradiance. Fourier series is also a good predictor of the 
cyclical dependence of solar radiation, combining several signifi- 
cant frequencies in its depiction. Boland [16,17] has shown that it 
captures yearly and intra-day cycles and can be used to effectively 
model the daily profile of solar irradiance time series, as well as 
capturing the variation over the year. 
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To judge the goodness of different detrending models, Ji and 
Chee [18] use the Augmented Dickey-Fuller (ADF) test to measure 
the stationarity of the detrended series. ADF is a test for unit root 
in a time series. 

A time series possesses a unit root if it is a realization of the 
process 


Xt =aXt1 +Z: (5) 


where a=1 and Z;~WN(0, 02). WN indicates white noise, that is 
independent and identically distributed. 

If there is a unit root in a time series, it is not stationary; There 
are other conditions for stationarity as well, for instance a cannot 
be greater than unity, but if all conditions hold then the time series 
may be considered stationary. 


2.1.3. ARMA model 

The ARMA model is based on two elementary models: the 

moving average (MA) model and the autoregressive (AR) model as 
p q 

SO= X aSt-) + X Aet- (6) 
i= j= 

In Eq. (6), S(t) is the forecasted solar irradiance at time t. In the 
AR model, p is the order of the AR process, and a; is the ith AR 
coefficient. In the MA model, q is the order of the MA error term, 2; 
is the jth MA coefficient and e(t) is the white noise that is 
uncorrelated random variables with zero mean and constant 
variance [19]. 

The Autoregressive Moving Average (ARMA) model is usually 
applied to auto correlated time series data. This model is a great tool 
for understanding and predicting the future value of a specified time 
series. ARMA is based on two parts: autoregressive (AR) part and 
moving average (MA) part. Also, this model is usually referred to as 
ARMA (p, q). In this p and q are the order of AR and MA respectively. 
The popularity of the ARMA model is its ability to extract useful 
statistical properties and the adoption of the well-known Box — 
Jenkins methodology (see Boland [17]). ARMA models are very 
flexible since they can represent several different types of time 
series by using different order. It has been proved to be competent in 
prediction when there is an underlying linear correlation structure 
lying in the time series. One major requirement for ARMA model is 
that the time series must be stationary (see Hamilton [20}). 


2.1.4. ARIMA techniques 

An extension of the ARMA models, the Auto-Regressive Inte- 
grated Moving Average (ARIMA) time series models form a general 
class of linear models that are widely used in modelling and 
forecasting time series in (see Box and Jenkins [19]). The ARIMA(p, 
d, q) model of the time series {X,,X2,...} is defined as 


@p(B)A'X; = Oq(B)ar (7) 
where 

®,(B) = 1—, B—B” —cldots—,B? (8) 
0,(B) = 1—0, B—02B? —---—0,B". (9) 


B is the backward shift operator, BXy=Xy_;, A=1-B is the 
backward difference, and , and @, are polynomials of order p and 
q, respectively. ARIMA(p, d, q) models are the product of an 
autoregressive part AR(p) (Eq. (8)), an integrating part I(d) = 47$, 
and a moving average MA(q) part (Eq. (9)). The parameters in ® 
and @ are chosen so that the zeros of both polynomials lie outside 
the unit circle in order to avoid generating unbounded processes. 

For some random shocks (being events such that the time series 
value departs from the time series average) drawn from a fixed 
distribution with mean zero and variance og, a sequence 
At, Q¢_1,t-2..., is called the white noise process. The backshift 


operator and the white noise process describe the intrinsic features 
of the time series: adjacent observations are dependent. a, is an 
independent shock of time step t of the white noise process. 

The ARIMA (Auto-Regressive Integrated Moving Average) tech- 
niques (see Hamilton [20]) are reference estimators in the predic- 
tion of global irradiance field. It is a stochastic process coupling 
autoregressive (AR) component to a moving average (MA) com- 
ponent, after differencing at appropriate time steps to remove any 
trends. It is in this way that ARIMA models allow treatment of 
non-stationary series. 

Reikard [5] applies a regression in log to the inputs of the 
ARIMA models to predict the solar irradiance. He compares ARIMA 
models with other forecast methods such as ANN. At the 24-h 
horizon, he states that the ARIMA model captures the sharp 
transitions in irradiance associated with the diurnal cycle more 
accurately than other methods. 


2.1.5. CARDS model 

Jing et al. [21] have developed a coupled autoregressive (AR) 
and dynamical system model to forecast solar radiation time series 
on hourly and intra-hourly time scales. 

The equation of the dynamical system part, the Lucheroni 
model is 


R=z (10) 


eZ = K(Z + R)—A(3R2Z + R?)—ez—yR—b + ¢ (11) 


where x,4,c,y and b are adjustable parameters, and ¢ is the 
noise term. 

In Eq. (10), R denotes the derivative of R with respect to time, 
and hence z stands for the second derivative of R with respect 
to time. 

The following discretized version of the model for our desea- 
soned solar radiation time series R; is used: 


Rey = Rt + 2 At + or (12) 


A 
Zt+1 =Zt + [ke + R)—åÀ (3R?z: + R?) —ezı—yRı—b] : = + dt (13) 


œt and a, are noise terms, A; is the time step. The parameters 
x,A,e,y and b can be estimated using the method of ordinary least 
squares (OLSs). 

The deseasoning was performed using Fourier series techni- 
ques as in Boland [16,17]. The residual series formed by subtract- 
ing the Fourier series component from the original series cannot 
be modeled sufficiently using and autoregressive process alone. 
This is because the autoregressive process is too efficient at mean 
reversion to be able to reach the peaks in the series. The 
introduction of a resonating model introduced for the power 
market by Lucheroni [22] plus the judicious intermittent use of a 
proxy for curvature allows for a much superior fit to this residual 
series. The results found for the CARDS model compare very 
favorably with what Kostylev and Pavlovski [4] found from their 
survey of the literature. From that paper, the best performing 
model at the 1h time step had a rRMSE of 17% for mostly clear 
days and 32% for mostly cloudy, whereas for the CARDS model it is 
16.5% for all days. 


2.2. Non-linear models 


For about one decade, there has been great interest in research 
on artificial intelligence (AI) techniques, not only for forecasting 
but also for a broad range of applications, including control, data 
compression, optimization, pattern recognition, and classification. 
An overview of the application of Al techniques for modeling and 
forecasting of the solar irradiance is presented in Mellit and 


M. Diagne et al. / Renewable and Sustainable Energy Reviews 27 (2013) 65-76 69 


Kalogirou [23] where several approaches have been compared and 
analyzed. 


2.2.1. Artificial neural network (ANN) 

As an alternative to conventional approaches, ANNs have been 
successfully applied for solar irradiance estimation. 

The most popular form of neural network is the so called 
multilayer perceptron (MLP) structure (see Lauret [24]). The MLP 
structure consists of an input layer, one or several hidden layers 
and an output layer. The input layer gathers the model's inputs 
vector x while the output layer yields the model's output vector y. 
Fig. 1 represents a one hidden layer MLP. 

The hidden layer is characterized by several non-linear units 
(or neurons). The non-linear function (also called activation 
function) is usually the tangent hyperbolic function f(x). 


ex—e* 


ex + e~ a9 


f= 


Therefore, a neural network with d inputs, h hidden neurons 
and a single linear output unit defines a non-linear parameterized 
mapping from an input x to an output y given by the following 
relationship: 


h d 
y=y(x; w= 2 ws 2,0) (15) 


j=0 


The parameters of the NN model are given by the so called 
weights and biases that connect the layers between them (notice 
that in Eq. (15), the biases are denoted by the subscripts i=0 and 
j=0 and are not represented in Fig. 1). The NN parameters, denoted 
by the parameter vector w, govern the non-linear mapping. 

ANNs recognize patterns in data and have been applied to solar 
forecasting. Using training data, ANNs reduce normalized root 
mean square error (rRMSE) of daily average GHI by as much as 15% 
when compared to 12-18 h ahead NWP forecasts (see Guarnieri 
et al. [25]). 

Heinemann et al. [7] use satellite images for horizons below 
6h; In Lorenz et al. [26] longer horizons of forecast produces by 
NWP models are used as input to an ANN to predict global 
irradiance. 

Mellit and Pavan [27] developed a Multilayer Perceptron MLP- 
model to forecast the solar irradiance 24h ahead. The proposed 
model accepts as input parameters mean daily irradiance and 
mean daily air temperature; The output is solar irradiance data 
24h ahead. Performance prediction of a grid-connected PV plant 
at Trieste, Italy, had a correlation coefficient of more than 98% for 
sunny days and slightly less than 95% for cloudy days. 


h hidden units 


Fig. 1. Sketch of a MLP with d inputs and h hidden units. Extracted from Lauret 
[24]. 


Kemmoku et al. [11] used a multistage ANN to predict GHI of 
the next day. The input data to the network are the average 
atmospheric pressure, predicted by another ANN and various 
weather data of the previous day. Irradiance forecast by the 
multi-stage and the single-stage neural networks are compared 
with measured irradiance. The results show that the mean bias 
error (MBE) reduces from about 30% (by the single-stage) to about 
20% (by the multi-stage). 

Sfetsos and Coonick [12] use ANN to make one-step predictions 
of hourly values of global irradiance and to compare them with 
linear time series models that work by predicting the clearness 
index. They introduced an approach for forecasting hourly solar 
irradiance using various artificial intelligence based techniques 
(ANN and ANFIS). They also investigated other meteorological 
variables such as temperature, wind speed, and pressure. 
A comparison between the various models in terms of prediction 
error and training time indicated that the network trained with 
the Levenberg-Marquardt (LM) algorithm network was as the 
optimum prediction model. 

Mihalakakou et al. [28] developed a total solar irradiance time 
series simulation model based on ANN and applied it in Athens. 
The Neural Logic Network was identified as the model with the 
least error. It incorporates Logic Rules that produced an rRMSE of 
4.9% lower than the persistence approach. 

Fatih et al. [29] developed a time delay neural network (TDNN) 
model from general feed forward neural network to obtain the 
relationship between the input and output position in time series. 
Conventional ANN provides their response to the weighted sum of 
the current inputs. For TDNN, it extends the sum to a finite 
number of past inputs. In this way, the output provided by a given 
layer depends on the output of the previous layer's computed 
values based on the temporal domain of input values. Because of 
the very similar structure of the TDNN and the general MLP, back- 
propagation with some modifications can be applied to train the 
TDNN. The strength of this algorithm is its ability to model 
nonlinear series. With TDNN, there is no need to specify a 
particular model form, since the model is adaptively formed based 
on the features presented by the data. This data driven algorithm 
is suitable for many time series where no theoretical model is 
available. 


2.2.2. Wavelet neural network 

Wavelet transform has time-frequency localization property 
and focal features and neural network (NN) has self-adaptive, fault 
tolerance, robustness, and strong inference ability. How can one 
combine the advantages of wavelet transform and NN to solve 
practical problems? The so-called wavelet neural network (WNN) 
or wavelet network (WN) is a variety of two techniques and 
inherits the advantages of the neural network and wavelet 
transformation. WNN uses the wavelet function as the activation 
function instead of the Sigmoid or tangent hyperbolic activation 
function. For WNN, the transfer function of hidden layer nodes is 
the mother wavelet function; and the network signal is prior to 
transmission while error is backpropagation in the training pro- 
cess. The network topology is shown in Fig. 2. 

In Fig. 2, X1,X2,...,Xn is the input vector; y;,y2,...,¥, is the 
predicted output; and w;j and wxj are the weights connecting every 
layer and hj is mother wavelet function. 

For the input signal sequence x = (xX1,X2,...,Xn), the output of 
the hidden layer is calculated as 


hG) = hj eer). j=1,2,...m, (16) 
J 


where h(j) is output value for the node j in the hidden layer; hj is 
the mother wavelet function; wy is weight connecting the input 
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Fig. 2. Topology of wavelet neural network. Extracted from Gaige [30]. 


layer and hidden layer; b; is the shift factor, and a; is the stretch 
factor for hj. 
The output of the output layer is calculated as 


m e 
y= X orh, 


t= 


k=1,2,...,1, (17) 


where h(i) is the output value for node i in the hidden layer; wj; is 
weight connecting the hidden layer and output layer; | and m are 
the number of nodes for output layer and the hidden layer, 
respectively. 

Mellit et al. [31] proposed an adaptive wavelet-network model 
for forecasting daily total solar irradiance. In this study, several 
structures have been investigated for resolving the missing data 
problem. In this particular estimation process, the model consists 
of an adaptive neural-network topology with the wavelet trans- 
formation embedded in the hidden units. 

Cao and Lin [32] proposed a new model for forecasting global 
solar irradiance based on diagonal recurrent wavelet neural net- 
work (DRWNN) and a special designed training algorithm. Simula- 
tion examples proved that the model is capable of mapping solar 
irradiance that is usually highly non-linear and time-changeable. 
This is because the DRWNN combines the advantages of both RNN 
(recurrent neural network) and wavelet neural network (WNN). 


2.2.3. ANN and classical time series models comparison 

A comparison of ANN and classical time series models has been 
carried out in Reikard [5] and in Sfetsos and Coonick [12]. Both 
studies find that the error of a simple regression model can be 
reduced considerably by a factor in the range of 0.6-0.8 when 
using advanced models. Reikard [5] compares a regression model, 
the UCM Model, an ARIMA, a transfer function model, a neural 
network model and hybrid model. For this study, the author uses a 
logarithmic scaling of the input. Results show that for the resolu- 
tions of 60, 30 and 15 min, the ARIMA model shows better results. 
In Sfetsos and Coonick [12], a feed-forward ANN is identified as the 
most appropriate. The analysis in Reikard [5] for several stations 
with different climatic conditions also shows that there is a strong 
influence of the climatic conditions on both forecast accuracy and 
potential for improvement by the use of advanced models. 


3. Cloud imagery and satellite based models 


The method from estimated solar radiation from satellite 
images is complex, involving several steps. For full details, please 
see Lorenz [3] and Rigollier [33]. 


3.1. Cloud imagery 


Besides the deterministic daily and annual patterns of irradi- 
ance, clouds cover as well as cloud optical depth have the 
strongest influence on solar irradiance at surface level. Clouds 
show a strong variability in time and space. Hence, determination 
of clouds at a designated time is an essential task in irradiance 
forecasting and modeling. For forecast horizons up to some hours, 
the temporal change of cloud structures is strongly influenced by 
cloud motion as a result of horizontal advection. 

Satellites and ground-based sky images have been used for the 
determination and forecasting of local solar irradiance conditions. 
The basis of this method relies upon the determination of the 
cloud structures during the previous recorded time steps. Extra- 
polation of their motion leads to a forecast of cloud positions and, 
as a consequence, to the local radiation situation. Satellites and 
ground-based sky images with their high temporal and spatial 
resolution offer the potential to derive the required information on 
cloud motion. 

Through processing of satellite or ground images, clouds can be 
detected, characterized, and advected to predict GHI relatively 
accurately up to 6h in advance. The time series models based on 
satellite data and sky images detect the motion of cloud structures 
using motion vector fields (see Lorenz et al. [6]). 

Irradiance for all sky conditions including cloudy skies may be 
derived using radiative transfer models (RTM) (see Heinemann 
et al. [7]) requiring input on the vertical structure of cloud physical 
parameters, for example, cloud and ice water content or droplet 
radius. Numerical weather prediction (NWP) models (see Section 
4) imply parameterizations of radiative transfer calculations. 

The errors of satellite data and sky images based forecasts 
proposed in the literature increase drastically under low sun 
elevations, high spatial variabilities and low irradiance conditions. 
Hammer et al. [34] demonstrated achieved 17% rRMSE in satellite 
imagery for 30 min cloud index forecasts and 30% rRMSE at 2h 
forecast horizons. For intra-day forecasts, a reduction in rRMSE by 
7% to 10% compared to persistence forecasts was found. 

Chow et al. [35] presented a technique for intra-hour, sub- 
kilometer cloud shadow nowcasting and forecasting using a ground- 
based sky images for selected days at the UC San Diego. This technique 
allows one to obtain sky cover, cloud motion, cloud shadows, 
irradiance, and to forecast cloud locations. 


3.2. Satellite images 


An approach to forecast solar irradiance based on Meteosat 
satellite images as a basis for PV power forecast was proposed in 
Lorenz et al. [36]. They investigated and compared various 
methods to derive motion vector fields from Meteosat data, and 
applied them to forecast solar irradiance up to some hours ahead. 
In Perez et al. [1], the results of irradiance forecasts based on the 
images of the Geostationary Operational Environmental Satellite 
(GOES) with a similar approach patterned after Lorenz et al. [6] 
are shown. 


3.3. Ground-based sky images 


To achieve high temporal and spatial resolution for intra-hour 
forecasts, NWP and satellite forecasts are currently inadequate. 
Ground observations using a total sky imager (TSI) present an 
opportunity to fill this forecasting gap and deliver a sub-kilometer 
view of cloud shadows over a large-scale PV power plant or an 
urban distribution feeder. 

Compared with satellite data, ground-based sky images offer 
a much higher spatial and temporal resolution, including the 
possibility of capturing sudden changes in the irradiance, often 


M. Diagne et al. / Renewable and Sustainable Energy Reviews 27 (2013) 65-76 71 


referred to as ramps, on a temporal scale of less than 1 min. The 
maximum possible forecast horizon strongly depends on the cloud 
speed and is limited by the time until the monitored cloud scene 
has passed the location or area of interest. This time is determined 
by the spatial extension of the monitored cloud scenes in combi- 
nation with cloud velocities. In Chow et al. [35], forecasts up to 
5 min ahead were evaluated for 4 partly cloudy days. An estima- 
tion of a maximum possible extension of the forecast horizon in 
dependence on the cloud scene resulted in values ranging from 
5 to 25 min. 

Only short deterministic forecast horizons are feasible using a 
single TSI at a site due to low clouds and large clouds variabilities 
at the fine spatial scale studied. Capturing these features determi- 
nistically is nearly impossible with satellite or NWP approaches 
(Chow et al. [35]). 


4. Numerical weather prediction models 


Numerical weather prediction (NWP) models are operationally 
used to forecast the state of the atmosphere up to 15 days ahead. 
The temporal development of the state of the atmosphere is 
modeled by the basic differential equations that describe the 
physical laws governing the weather [3]. 

Starting from initial conditions that are derived from world- 
wide observations, in a first step, the future state of the atmo- 
sphere is calculated with a global NWP model. Global NWP models 
are currently in operation at about 15 weather services. Examples 
are the Global Forecast System (GFS) run by the US National 
Oceanic and Atmospheric Administration (NOAA) and the 
European Centre for Medium-Range Weather Forecasts (ECMWF). 
Global models usually have a coarse resolution and do not allow 
for a detailed mapping of small-scale features, although resolution 
has increased rapidly during the last few years and nowadays, 
depending on the model, is in the range of 16-50 km. In the next 
step, different concepts may be applied to account for local effects 
and to derive improved site-specific forecasts. One possibility is 
the downscaling by mesoscale models, which are also referred to 
as regional models. Mesoscale models cover only a part of the 
Earth but can be operated with a higher spatial resolution. They 
are routinely run by national weather services and private weather 
companies. Also, postprocessing methods, may be applied to 
model local effects. They allow the correction of systematic 
deviations in dependence on different meteorological parameters 
and for modeling of the irradiance if it is not provided as output 
parameter of an NWP model. In the next few subsections we will 
review NWP configuration, input data feature, global model 
example, mesoscale model example, NWP accuracy, NWP limita- 
tion and postprocessing methods. 


4.1. NWP configuration 


Before running an NWP model, it should be configured. The 
principal variables to fix are the time step of internal calculation 
and the horizontal resolution of grid points. 


4.1.1. Temporal resolution 

The internal time step gives the period over which the change 
of the atmospheric variables is described by the dynamic equa- 
tions. Temporal resolution of internal calculations in NWP models 
usually is considerably higher than that of the output variables. 
Output variables are delivered with a resolution of typically 1 h for 
regional models and 3-6 h for global models. This internal time 
step may be down to 30s for highly resolved calculations with 
mesoscale models and is about 10 min for global NWP models. 


4.1.2. Spatial resolution 

The horizontal resolution determines the spatial extent of 
weather phenomena that can be directly simulated. Grid points 
are usually distributed equally in the horizontal range. The resolu- 
tion of global NWP models nowadays is in the range of 16-50 km. In 
mesoscale models, the horizontal resolution may be down to 1 km; 
weather services typically operate mesoscale models with a spatial 
resolution in the range of 5-20 km. The resolution of vertical levels 
is generally adapted to the occurrence of physical processes that take 
place in certain regions of the atmosphere. 


4.2. Input data features 


To start a forecast, information on the current state of the 
atmosphere is necessary. For global NWP models, this information 
is obtained from a worldwide network of meteorological observa- 
tions and measurements. The key variables needed are the three- 
dimensional fields of wind, temperature, and humidity and the 
two-dimensional field of surface pressure. Boundary variables like 
snow cover or sea surface temperature are also of high impor- 
tance. Regional models use initial conditions as well as lateral 
boundary conditions from global NWP model output, and also 
offer the possibility of integrating local measurements. 


4.3. Global model example: ECWMF 


ECMWF provides weather forecasts up to 15 days ahead, 
including solar surface irradiance and different cloud parameters 
as model output. ECMWF forecasts have shown their high quality 
as a basis for both wind and solar power forecasts. These forecasts 
are described here as an example of global NWP model forecasts. 

The evaluations of ECMWF-based irradiance in Lorenz et al. 
[36-38] are based on the T799 version with a spatial resolution of 
25 km x 25 km. The current version T1279 was implemented in 
January 2010 and shows a horizontal resolution of 16 km x 16 km. 
Ninety-one hybrid vertical levels resolve the atmosphere up to 
0.01 hPa corresponding to approximately 80 km. The temporal 
resolution of the forecasts is 3 h for the first 3 forecast days that 
are most relevant for PV power prediction. 


4.4, Mesoscale models example: MM5 and WRF 


The non-hydrostatic, fifth-generation mesoscale model MM5 
has been developed at Pennsylvania State University and subse- 
quently at the National Center for Atmospheric Research (NCAR). It 
uses a terrain-following coordinate, solves its finite-difference 
equations with a time-split scheme, and has multiple nesting 
capabilities (Grell et al. [39]). The WRF model is designed to be a 
flexible, state-of-the-art model and is developed as a collaborative 
effort of several institutes. WRF is supported as a community 
model with continuous development and integrates features of 
different mesoscale models, including also MM5 and the Eta 
model of the National Centers for Environmental Prediction 
(NCEP). In this sense, WRF can be seen as a follow-up model 
to MM5. The current WRF model, version 3, is described in 
Skamarock et al. [40]. Both WRF and MM5 offer a number of 
parameterizations for the different physical processes. This allows 
adapting the configuration of the model to the specific climatic 
conditions for an interested region. In addition, the capability of 
MM5 and WRF to integrate local measurements, for example, 
aerosols, may also contribute to improving forecast accuracy. The 
simulation of mesoscale and small-scale phenomena, which is 
essential for calculations with high spatial resolution, is supported 
by the non-hydrostatic dynamics. With the possibility of high 
spatial resolution, the effects of topography may be considered in 
much more detail than for large-scale models. 
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Fig. 3. Input data to WRF model (Skamarock et al. [40]). 


4.4.1. Input data 

Mesoscale models require input from global NWP models for 
initialization and boundary conditions. Frequently, GFS data of NOAA 
are used to initialize MM5 or WRF for operational applications, 
because, in contrast to ECMWF data, they are available for free. The 
input data used have a significant influence on the results, especially 
for cloudy conditions. This has been demonstrated in a case study 
reported in Heinemann et al. [7], where different NWP models 
(ECMWF and the global and local models of the German Weather 
Service (DWD)) were compared for initialization. Fig. 3 illustrates 
patterns of input data for regional NWP model. 


4.4.2. Grid nesting 

To achieve the intended high spatial resolution in a mesoscale 
model with reasonable computing time, the resolution of the 
driving global model is increased stepwise with internal nesting. 
For example, in the study reported in Heinemann et al. [7], the 
outer domain of MM5 covering large parts of Europe has a 
resolution of 27 km x 27 km, the next domain has a resolution of 
9 km x 9 km, and the final resolution of the innermost domain is 
3 km x 3 km. 


4.4.3. MM5 and WRF configuration 

Recently, several research groups have investigated the potential of 
MM5 and WRF irradiance forecasts for solar energy applications. In a 
first step, an appropriate setup of MM5 or WRF has to be determined. 
Heinemann et al. [7] has addressed this task by comparing different 
configurations of the cumulus, moisture, and planetary boundary layer 
parameterizations with respect to irradiance calculations in a case 
study. To limit the computational effort when rerunning simulations 
for several configurations, a set of 6 test days was defined covering 
different cloud conditions: clear sky, broken clouds, and overcast. 

Heinemann et al. [7] show an evaluation of MM5 forecasts for a 
40-day period in summer 2003 in Southern Germany in comparison 
with other forecasting approaches. Two studies comparing different 
methods to predict solar irradiance including WRF forecasts for 
different locations in the United States are reported in Remund et al. 
[41] and in Perez et al. [42]. In Lara-Fanego et al. [43], a detailed 
evaluation study of WRF irradiance forecasts in Andalusia (Southern 
Spain) is given with calculations of 1 month for each season. 


4.5. NWP model accuracy 


Forecasts beyond 6h, up to several days ahead, are generally 
most accurate if derived from NWP models. NWP models predict 
GHI using columnar (1D) radiative transfer models (RTM) (see 
Heinemann et al. [7]). Heinemann et al. [7] showed that the MM5 
mesoscale model can predict GHI in clear skies without bias. 
However, the bias was highly dependent on cloud conditions 
and becomes strong in overcast conditions. 

Perez et al. [44] examined the accuracy of the National Digital 
Forecast Database (NDFD), a derivative of the operational NWP 


models published by the NCEP. After a local correction function was 
applied, results show that for 8-26 h forecast horizons, the NDFD 
had an hourly average GHI relative RMSE (rRMSE) of 38%. 

Remund et al. [41] evaluated different NWP-based GHI forecasts in 
USA, reporting rRMSE values ranging from 20% to 40% for a 24h 
forecast horizon. Similar results were reported by Perez et al. [1], 
evaluating NWP-based irradiance forecasts in several places in the 
USA. Remund et al. [41] examined NWP biases compared to a single 
site and find that ECMWF and GFS next day GHI forecasts have a Mean 
Bias Error (MBE) of 19%. This MBE was found to be approximately 
constant for intra-day (hour-ahead) to 3 days ahead forecast horizons. 

Lorenz et al. [45] evaluated several NWP-based GHI forecasts in 
Europe. Overall, results showed rRMSE values of about 40% for 
Central Europe and about 30% for Spain. Evaluating ECMWF accu- 
racy in Germany, Lorenz et al. [36] showed that NWP MBE was 
largest for cloudy conditions with moderate clear sky index 
(0.3 < kt* < 0.6), while forecasted clear conditions were relatively 
unbiased. They reported rRMSE values of about 35% for single 
stations for a 24h horizon forecasts. 


4.6. NWP limitations 


A limitation of NWP forecasting is its coarse resolution. Even the 
0.1° x 0.1° NAM spatial resolution is insufficient to resolve most clouds. 
Only an average cloud cover can be forecasted for a given point. For 
global models (GFS and ECMWF) the resolution is even coarser. 
However, even if the spatial resolution is finer, the temporal output 
intervals would not permit the assessment of time dependent cloud 
cover variability, important in predicting ramp rates and ranges of 
variability for solar power plants. Although NWP model time-steps are 
on the order of minutes, the RTM are run less frequently, and the 
output is only hourly (NAM) or every 3h (GFS and ECMWF). 
Consequently, any patterns with characteristic time scales less than 
an hour are unresolved. Linking observed temporal variability in GHI 
to native NWP forecasts will require further research. 


4.7. Postprocessing methods 


Postprocessing methods are frequently applied to refine the output 
of NWP models like presented in Lorenz [3]. In particular, they may be 
utilized to: 


® reduce systematic forecast errors (correction of systematic 
deviations); 

© account for local effects (e.g., topography); 

è derive parameters that are not directly provided by the NWP 
models (e.g., solar surface irradiance is still not a standard 
output parameter); 

© combine the output of different models in an optimum way. 


Various approaches have been proposed to address these 
issues, some of which are presented in the following subsections. 
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4.7.1. Model output statistics (MOS) 

MOS relates observed weather elements to appropriate vari- 
ables (predictors) via a statistical approach. These predictors may 
be NWP model forecast, prior observations, or geoclimatic data. 

A state-of-the-art MOS for solar irradiance predictions based on 
ECMWF forecasts has been introduced in Bofinger and Heilscher 
[46]. Multiple regression is applied to modify long-term monthly 
mean values of the forecasted data. Direct model output of 
ECMWF and statistically derived predictors are used to create 
daily solar electricity predictions accurate to 24.5% rRMSE for 
averaged daily forecasts. The MOS is operated on the basis of 
ground-measured irradiance values when available. For locations 
without irradiance measurements, irradiance derived from Meteo- 
sat data with the Heliosat method is used instead. A comparison of 
irradiance forecasts using this MOS scheme with WRF forecasts 
and other approaches is given in Heinemann et al. [7] and Lorenz 
et al. [36]. 

Lorenz et al. [36] related forecasted solar zenith angle (SZA) and 
clear sky index to ECMWF MBE for Germany, revealing a consis- 
tent over-prediction (up to 100 W m~?) for moderately cloudy 
conditions. Using a MOS correction function eliminated bias and 
reduced RMSE of hourly forecasts by 5% for 24h forecasts. A 
stepwise multivariate fourth-order regression was applied to 
derive the MOS correction function. In Mathiesen and Kleissl 
[47] the analysis and MOS correction of GHI forecasts from NAM, 
GFS, and ECMWF models within the continental United States was 
presented. They indicated that MOS application to the NWP 
irradiance output was successful in minimizing bias and reducing 
RMSE, but did not provide information about bias source. MOS 
corrections in the measured clear sky regime did not reduce RMSE. 
This is because the MOS could not distinguish between RTM errors 
(over-prediction of GHI even for clear skies, especially for NAM) 
and cloud model errors (incorrect parameterization of RTM 
inputs). Consequently, many initially accurate forecasts were 
unnecessarily corrected. Differentiating between the sources of 
the error is important to selectively correct forecasts. Although 
traditionally MOS schemes are mostly based on linear regression, 
any statistical approach relating observed variables to NWP output 
fits to the concept of MOS. In particular, ANN have also been used 
to improve NWP output with respect to irradiance prediction (Cao 
and Lin [48], Guarnieri et al. [25]). In Guarnieri et al. [25], an ANN 
is applied to irradiance forecasts of the NCEP Eta model run 
operationally at the Brazilian Center for Weather Forecasts and 
Climate Studies (CPTEC/INPE). An evaluation with measurements 
for two stations in the south of Brazil reveals a strong over- 
estimation of the irradiance by the original forecasts, and a 
considerable improvement is achieved by the application of an 
ANN using different atmospheric forecast parameters of the Eta 
model as input. 


4.7.2. Kalman filter 

Systematic deviations of NWP output variables often depend 
on the meteorological situation. In Lorenz et al. [36], a bias 
correction in dependence on the predicted cloud situation for 
the application to ECMWF irradiance forecasts is introduced. 
A method for bias removal of irradiance forecasts using Kalman 
filtering is introduced and compared with the bias correction 
according to Lorenz et al. [36]. Kalman filters are designed to 
efficiently extract a signal from noisy data and are therefore 
expected to show a more robust performance if only limited 
training data are available, which is the case if the training is 
performed on the basis of individual stations. 

Pelland et al. [49] found that the most suitable realization of 
their approach was a set of Kalman filter equations established 
separately for each forecast horizon and modeling the bias in 


dependence on the forecasted irradiance. The accuracy assessment 
was performed for single stations and for regional average values. 
At the level of individual stations, the bias removal based on 
Kalman filtering outperforms the other approach. However, the 
improvement compared with the original forecasts is small for 
single stations, while for regional averages both bias removal 
approaches significantly reduced the RMSE. 


4.7.3. Temporal interpolation 

Global model forecasts are provided with a temporal resolution 
of 3-6h. The management of electricity grids, however, needs 
forecasts of the expected solar power input at least on an hourly 
basis. Different interpolation techniques may be applied in order 
to derive hourly forecasts from global NWP output. Lorenz et al. 
[36] propose an approach combining the forecast data with a 
clear-sky model to account for the typical diurnal course of 
irradiance. 


4.7.4. Spatial averaging 

Lorenz et al. [36] found that for ECMWF forecasts with an 
original spatial resolution of 25 km x 25 km and temporal resolu- 
tion of 3h, best results are achieved for average values of 4 x 4 
grid points corresponding to a region of 100 km x 100 km (Lorenz 
et al. [36]). However, the improvement compared with forecasts 
that evaluate only the next grid point is small due to the already 
coarse spatial and temporal resolution of the original forecasts. 

Spatial averaging has a much stronger impact for mesoscale or 
multi-scale model output with hourly values and a finer grid 
resolution. An analysis of the high-resolution WRF forecasts 
presented in Heinemann et al. [7] showed that averaging irradi- 
ance predictions over an area of 180 km x 180 km reduces the 
RMSE to approximately 85% of the RMSE when evaluating the 
nearest grid point only. Similar improvements are achieved for 
WRF forecasts provided by Meteotest that are delivered as average 
values of 10x10 model pixels, corresponding to an area of 
50 km x 50 km (Lorenz et al. [36]), and also for irradiance forecasts 
with the Canadian GEM model, where averaging areas in the range 
of 300 km x300km to 600kmx600km gave better results 
(Pelland et al. [49]). Mathiesen and Kleissl [47] report 100 km x 
100 km as a suitable averaging area for irradiance forecasts of GFS 
model and NAM model. 


4.7.5. Physical postprocessing approaches 

A few studies also investigate physical postprocessing proce- 
dures involving radiation transfer calculations. This allows for 
integrating additional parameters that are generally not modeled 
in detail with NWP models, for example, aerosols. A partly physical 
postprocessing procedure for topographic downscaling of solar 
irradiance forecasts in mountainous regions is proposed in 
Lara-Fanego et al. [43]. The disaggregation is carried out by 
accommodating the initial WRF irradiance estimates to the eleva- 
tion of a target digital elevation model with a spatial resolution of 
90 m x 90m. The proposed method accounts for shading, sky- 
view reduction, reflected irradiance, and scaling to the inclined 
terrain surface. 


4.8. Human interpretation of NWP output 


Finally, a traditional method to obtain improved local forecasts 
from NWP model output is the participation of a human forecaster 
[16]. Meteorologists at weather services routinely analyze and com- 
pare the output of different global and local NWP models and 
meteorological measurements. In particular, they also use their expert 
knowledge to decide on the final forecast values, for example, of cloud 
cover. Solar irradiance forecasts may be derived by combining the 
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cloud cover forecasts of meteorologists with a clear-sky model (see 
Section 2.1.2). An advantage of this approach is that forecasts may be 
adjusted for local events or weather situations difficult to forecast with 
NWP models or statistical methods, like, fog. 


5. Hybrid models 


Hybrid models have been introduced to overcome the deficiency 
of using a individual model such as statistical methods (ARIMA, 
Multiple Regression, etc.) and AI methods. Hybrid models merge 
different methods to improve the prediction accuracy. Hybrid models 
can be also referred as combined models or ensemble models and 
often these terms are used synonymously. Hybrid methods can be 
implemented in three different ways; linear models, nonlinear 
models and both linear and nonlinear models. 

With the intention to improve the forecasting accuracy, the 
combination of forecasting approaches has been proposed by many 
researchers [50,13,32,18,5]. From their studies, they indicate that the 
integrated forecasting techniques outperform the individual forecasts. 

Artificial intelligence techniques, such as fuzzy logic and neural 
networks, have been used for estimating hourly global solar irra- 
diance from satellite images. The results seem to point out that fuzzy 
lo>gic and neural network models are better than regression models. 

Cao and Cao [50,13] developed a hybrid model for forecasting 
sequences of total daily solar irradiance, which combines ANN with 
wavelet analysis. 

Cao and Lin [32] use an ANN (with a special designed training 
algorithm) combined with wavelets (based on diagonal recurrent 
wavelet neural network (DRWNN)) to predict next day hourly values 
of global irradiance. Different types of meteorological observations 
are used as input to the models; among others the daily mean global 
irradiance and daily mean cloud cover of the day are forecasted. 

Ji and Chee [18] use a hybrid model of ARMA and TDNN to 
improve the prediction accuracy. They suppose that the daily solar 
irradiance series is composed by linear and nonlinear components 
and use the ARMA model to fit the linear component and the 
TDNN model to find the nonlinear pattern lying in the residual. 
This hybrid model has the potential to harness the unique features 
and strengths of both models. It is more accurate than using the 
ARMA or TDNN models separately. 
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6. Future solar irradiance forecasting approaches 
for small-scale insular grids 


A summary of the literature on solar irradiance forecasting models 
illustrated by Figs. 4 and 5 gives indications for future work. Fig. 4 
shows the classification of forecasting models based on spatial 
resolution of input data and temporal resolution of output or foreseen 
data. Fig. 5 illustrates the relation between the forecasting horizons, 
the forecasting models and the related activities for grid operators. 

However, a further consideration in choosing among forecast- 
ing models is efficiency. Reunion Island is a small territory with a 
high relief, a lot of microclimates and cloud formation processes. 
In this context, the selection of forecasting model is based on small 
horizontal and temporal resolution of forecasting model. 

Based on these consideration, for day ahead and intra-days 
forecast horizon, we find that global models like ECMWF and GFS 
whom present generally reliable results are limited by their coarse 
resolution for Reunion Island. We suggest to use the mesoscale 
model WRF. It allows for small horizontal and temporal resolution. 

For intra-day and intra-hour forecast horizon time series models 
are used. The large number of data (GHI) measured on the ground 
offers a large set of temporal series of irradiance. This time series will 
permit building a statistical forecasting model. In this forecast 
horizon, the ARIMA models seem to be the most reliable model. 
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Fig. 5. Relation between horizons, models and activities. 
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Fig. 4. Classification of model based on spatial and temporal resolution. 
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They can provide a forecast in a fraction of a second on a personal 
computer. However, for the horizon of forecast of a few minutes, the 
persistence model achieves better accuracy than ARIMA models. In 
this sense, the choice of the model depends critically on the horizon 
of forecast. At longer horizons, the data are dominated by the diurnal 
cycle. In this case ARIMA models work better. At higher frequency, 
the data is more dominated by short-term patterns which can be 
picked up by persistence or ANN. 


7. Conclusions 


Solar irradiance forecasting is important for the integration of 
photovoltaic plants into an electrical grid. Proper solar irradiance 
forecasting helps the grid operators to optimize their electricity 
production and /or to reduce additional costs by preparing an 
appropriate strategy. 

A number of time series models and numerical weather prediction 
(NWP) models have been reviewed in this paper. From the description 
of the various results of solar irradiance forecasting, we maintain that 
the choice of the appropriate forecasting models depends on forecast 
horizon and the available data. For forecast horizon from 6h up to 
3 days ECMWF associated with a MOS post-process shows the most 
accurate results. However, in the case of Reunion Island, the WRF 
model seems to be more pertinent. For a smaller forecast horizons, 
from 5 min to 4h ARIMA seems to present the best accuracy. Cloud 
imagery and a hybrid model can improve the results of forecasting 
when solar irradiance presents a strong variability like in many of 
insular territories. 

It is worth noting that Kostylev and Pavlovski [4] have done 
extensive analysis of the best performing models on differing time 
scales. As stated above, different models are best at different 
forecast horizons. They give the estimates of the best rRMSE 
values for both mostly cloudy times and mostly clear. The values 
range from approximately 17% at one hour, to 22% for three days 
for mostly clear. For mostly cloudy, they range from 33% to 44% 
over the forecast range of one to three days. 

Future work will include several elements to improve forecast 
accuracy. Sky image techniques will be used to account for the process 
of cloud formation. The interesting methods identified here (WRF, 
ARIMA and AR) will be combined to sky images to yield a compre- 
hensive and more accurate forecast product with different horizons of 
forecast. A post processing method such as Kalman filtering will be 
applied in order to improve accuracy of WRF forecasts. 

The method we propose will be faced with two limitations. The 
first is the need to have a high temporal and spatial resolution of the 
input data. The second, regarding the forecasts to be produced by the 
NWP models, is the need of high capacity computing power. Con- 
cerning post-processing methods, they are exclusively based on 
statistical properties of the time series and do not use the physical 
properties of solar irradiance. In the case of ARIMA models one 
limitation will be then to take into account of physical behavior of 
time series objects like sunrise and sunset. 

The goal is to take care of the needs of grid operators. 
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